/*
Replication Code for Figures 1 and 2 in the Text of "Stakeholder Participation 
in Policymaking: \\ Evidence from Medicare Fee Schedule Revisions"
*/
clear

cd "/Users/stevenrashin/Dropbox/Medicare_SG/Replication/"

set more off

use "figures_1_2.dta"

/******1. Binscatter analysis of commenting and stakes ******/
/*Note: We do this by hand rather than use binscatter so we can adopt a more 
flexible functional form than linear or quadratic fit, and also so we don't extrapolate where there's no data*/

/*Aggregate over all years*/
capture drop residy-n2015
capture drop n
gen residy=.
gen residx=.

capture drop residymass_spec-residxindiv_nospec n_all n
gen residymass_spec=.
gen residxmass_spec=.

/* Begin Figure 1 Code */

/* Figure 1 No Fixed Effects (Left Panel) */

/* First run the regressions then take the residuals for y*/

capture drop residyall -n
gen residyall=.
gen residxall=.

qui reg anycomment ln_total_allowed_sq hhi ln_medscape_salary_geo y_1 y_2 y_3
qui predict ry if e(sample), resid
qui summ anycomment if e(sample), meanonly
qui replace ry = r(mean) + ry
qui replace residyall = ry
drop ry

/* First run the regressions then take the residuals for x*/

qui reg stakes ln_total_allowed_sq hhi ln_medscape_salary_geo y_1 y_2 y_3
qui predict rx if e(sample), resid
qui summ stakes if e(sample), meanonly
qui replace rx = r(mean)+rx 
qui replace residxall = rx
drop rx

fastxtile n_all = residxall, nq(50)
gen n = .
replace n = n_all
preserve
collapse (mean) residyall (median) residxall stakes, by(n)
twoway 	   (scatter residyall residxall, mcolor(black) msymbol(circle_hollow)) ///
	       (lowess residyall residxall if residxall <= 0, lcolor(black) lpattern(dash)) ///
	       (lowess residyall residxall if residxall >=0, lcolor(black) lpattern(dash)) ///
		    ,  name(g1, replace) ///
			ytitle("Frequency (Residualized)") xtitle("Stakes (Residualized)")  title("No Fixed Effects") ///
			xlabel(-5 "-100K" -4 "-10K" -3 "-1K" -2 "-100" -1 "-10" 0 "0" 1 "10" 2 "100" 3 "1K" 4 "10K") ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) ylabel(0(0.005)0.015) scheme(plotplain)
restore
drop residyall -n

/* Figure 1 Fixed Effects (Right Panel) */
/* Same sequence as above */

capture drop residyall -n
gen residyall=.
gen residxall=.

local apolitical_indvars = "positive_stakes negative_stakes stake_dummy ln_total_allowed_sq hhi ln_medscape_salary_geo" 
local yeardummies = "y_1 y_2 y_3"
xtlogit anycomment `apolitical_indvars' `yeardummies', fe i(npi)
gen xt_touse = e(sample)

qui areg anycomment ln_total_allowed_sq hhi ln_medscape_salary_geo y_1 y_2 y_3 if xt_touse, absorb(npi)
qui predict ry if e(sample), resid
qui summ anycomment if e(sample), meanonly
qui replace ry = r(mean) + ry
qui replace residyall = ry
drop ry

qui areg stakes ln_total_allowed_sq hhi ln_medscape_salary_geo y_1 y_2 y_3 if xt_touse, absorb(npi)
qui predict rx if e(sample), resid
qui summ stakes if e(sample), meanonly
qui replace rx = r(mean)+rx 
qui replace residxall = rx
drop rx

fastxtile n_all = residxall, nq(50)
gen n = .
replace n = n_all
preserve
collapse (mean) residyall (median) residxall stakes, by(n)
twoway 	   (scatter residyall residxall, mcolor(black) msymbol(circle_hollow)) ///
	       (lowess residyall residxall if residxall <= 0, lcolor(black) lpattern(dash)) ///
	       (lowess residyall residxall if residxall >=0, lcolor(black) lpattern(dash)) ///
		    ,  name(g2, replace) ///
			ytitle("Frequency (Residualized)") xtitle("Stakes (Residualized)")  title("Physician Fixed Effects") ///
			xlabel(-5 "-100K" -4 "-10K" -3 "-1K" -2 "-100" -1 "-10" 0 "0" 1 "10" 2 "100" 3 "1K" 4 "10K") ///
			xline(0, lcolor(black) lpattern(dash)) legend(off)  scheme(plotplain)
restore
drop residyall -n

gr combine g1 g2, graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white))  xsize(5) ysize(2.5)

/* Begin Figure 2 Code */

/*Do for Mass Comments (Right Panel) */
drop residy residx residymass_spec residxmass_spe

capture drop residy-n2015
capture drop n
gen residy=.
gen residx=.

capture drop residymass_spec-residxindiv_nospec n_all n
gen residymass_spec=.
gen residxmass_spec=.

qui reg mass_comment ln_total_allowed_sq hhi ln_medscape_salary_geo y_* i.provider_id
qui predict ry if e(sample), resid
qui summ mass_comment if e(sample), meanonly
qui replace ry = r(mean) + ry
qui replace residymass_spec = ry
drop ry

qui reg stakes ln_total_allowed_sq hhi ln_medscape_salary_geo y_* i.provider_id
qui predict rx if e(sample), resid
qui summ stakes if e(sample), meanonly
qui replace rx = r(mean)+rx 
qui replace residxmass_spec = rx
drop rx

fastxtile n_spec = residxmass_spec, nq(50)
gen n = .
replace n = n_spec
preserve
collapse (mean) residymass_spec (median) residxmass_spec stakes, by(n spec_comment)
twoway 	   (scatter residymass_spec residxmass_spec if spec_comment == 1, mcolor(black) msymbol(circle)) ///
	       (scatter residymass_spec residxmass_spec  if spec_comment == 0, mcolor(black) msymbol(circle_hollow)) ///
		   (lowess residymass_spec residxmass_spec if residxmass_spec <= 0 & spec_comment == 1, lcolor(black) lpattern(solid)) ///
	       (lowess residymass_spec residxmass_spec if residxmass_spec >=0 &  spec_comment == 1, lcolor(black) lpattern(solid)) ///
	       (lowess residymass_spec residxmass_spec if residxmass_spec <= 0  & spec_comment == 0, lcolor(black) lpattern(dash)) ///
	       (lowess residymass_spec residxmass_spec if residxmass_spec >=0 &  spec_comment == 0, lcolor(black) lpattern(dash)) ///
		,  name(g3, replace) ///
			ytitle("Frequency (Residualized)") xtitle("Stakes (Residualized)")  title("Mass Comments") ///
			xlabel(-5 "-100K" -4 "-10K" -3 "-1K" -2 "-100" -1 "-10" 0 "0" 1 "10" 2 "100" 3 "1K" 4 "10K") ///
			xline(0, lcolor(black) lpattern(dash)) ///
			legend(order(1 2 4 6) lab(1 "Specialty Society Commented") lab(2 "Specialty Society Abstained") lab(4 "") lab(6 "") rowgap(0.1) ring(0) position(2) col(2) region(lstyle(solid))) ///
			ylabel(0(0.005)0.017) scheme(plotplain) xsize(2.5) ysize(2)
restore

/*Do for Tailored Comments (Left Panel)*/

drop residymass_spec-n


gen residyindiv_spec=.
gen residxindiv_spec=.


qui reg non_mass_comment ln_total_allowed_sq hhi ln_medscape_salary_geo y_* i.provider_id
qui predict ry if e(sample), resid
qui summ non_mass_comment if e(sample), meanonly
qui replace ry = r(mean) + ry
qui replace residyindiv_spec = ry
drop ry

qui reg stakes ln_total_allowed_sq hhi ln_medscape_salary_geo y_* i.provider_id 
qui predict rx if e(sample), resid
qui summ stakes if e(sample), meanonly
qui replace rx = r(mean)+rx 
qui replace residxindiv_spec = rx
drop rx


fastxtile n_spec = residxindiv_spec, nq(50)
gen n = .
replace n = n_spec
preserve
collapse (mean) residyindiv_spec  (median) residxindiv_spec  stakes, by(n spec_comment)
twoway 	   (scatter residyindiv_spec residxindiv_spec if spec_comment == 1, mcolor(black) msymbol(circle)) ///
	       (scatter residyindiv_spec residxindiv_spec if spec_comment == 0, mcolor(black) msymbol(circle_hollow)) ///
		   (lowess residyindiv_spec residxindiv_spec if residxindiv_spec <= 0 & spec_comment == 1, lcolor(black) lpattern(solid)) ///
	       (lowess residyindiv_spec residxindiv_spec if residxindiv_spec >=0 & spec_comment == 1, lcolor(black) lpattern(solid)) ///
	       (lowess residyindiv_spec residxindiv_spec if residxindiv_spec <= 0 & spec_comment == 0, lcolor(black) lpattern(dash)) ///
	       (lowess residyindiv_spec residxindiv_spec if residxindiv_spec >=0 & spec_comment == 0, lcolor(black) lpattern(dash)) ///
		,  name(g4, replace) ///
			ytitle("Frequency (Residualized)") xtitle("Stakes (Residualized)")  title("Tailored Comments") ///
			xlabel(-5 "-100K" -4 "-10K" -3 "-1K" -2 "-100" -1 "-10" 0 "0" 1 "10" 2 "100" 3 "1K" 4 "10K") ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) ylabel(0(0.005)0.017) scheme(plotplain) xsize(2.5) ysize(2)
restore
drop residyindiv_spec-n


grc1leg g4 g3, legendfrom(g3) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) name(gr_massspec, replace)
graph display gr_massspec, ysize(2.5) xsize(5)
